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Ab initio numerical simulations of the thermodynamics of dense quark matter remain a challenge. 

Apart from the infamous sign problem, lattice methods have to deal with finite volume and dis¬ 
cretization effects as well as with the necessity to introduce sources for symmetry-breaking order 
parameters. We study these artifacts in the Polyakov-loop-extended Nambu-Jona-Lasinio model, 
and compare its predictions to existing lattice data for cold and dense two-color matter with two 
flavors of Wilson quarks. To achieve even qualitative agreement with lattice data requires the in¬ 
troduction of two novel elements in the model: (i) explicit chiral symmetry breaking in the effective 
contact four-fermion interaction, referred to as the chiral twist, and (ii) renormalization of the 
Polyakov loop. The feedback of the dense medium to the gauge sector is modeled by a chemical- 
potential-dependent scale in the Polyakov-loop potential. In contrast to previously used analytical 
ansatze, we determine its dependence on the chemical potential from lattice data for the expectation 
value of the Polyakov loop. Finally, we propose to add a two-derivative operator to our effective 
model. This term acts as an additional source of explicit chiral symmetry breaking, mimicking an 
analogous term in the lattice Wilson action. 

I. INTRODUCTION 

Hadronic matter in extreme conditions such as high 
temperature or high density has received considerable at¬ 
tention over the past decades. However, direct numerical 
simulations of the theory of strong interactions — quan¬ 
tum chromodynamics (QCD) — at nonzero baryon den¬ 
sity are a formidable challenge due to the infamous sign 
problem. Large efforts have been made to overcome this 
problem (see for instance Ref. [1] for some recent propos¬ 
als), yet with limited success so far. At present, the only 
tools for quantitative analysis of dense nuclear matter are 
phenomenological effective models and, to some extent, 
continuum functional methods [2]. 

Lattice simulations and phenomenological models can 
be of mutual benefit: while numerical simulations can 
provide a hrm model-independent basis for effective con¬ 
tinuum approaches, results obtained by the latter can 
be easily extrapolated to conditions where lattice tech¬ 
niques are difficult to apply. However, there is a gap 
that needs to be bridged to make this interaction possi¬ 
ble. Lattice simulations have to deal with several arti¬ 
facts, namely the effects of finite volume and spacetime 
discretization, as well as the need for external sources to 
pick a unique ground state whenever continuous symme¬ 
tries are expected to be spontaneously broken. Usually, it 
is much easier to introduce these effects in models, rather 
than eliminate them from lattice simulations. 

The main goal of the present paper is to do exactly 
that. We effectively re-introduce the effects of exter¬ 
nal sources and spacetime discretization using a phe- 
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nomenological continuum model and focus on discrimi¬ 
nating between physics and lattice artifacts. Our work 
is based upon results of recent simulations of two-color 
QCD (2cQCD) at high baryon density, using two fla¬ 
vors of Wilson-type quarks [3, 4]. These are compared 
to an effective model of the Nambu-Jona-Lasinio (NJL) 
type [5-7], augmented with the Polyakov loop, which is 
an (approximate) order parameter for deconfinement [8- 
10]. We discuss in detail the model under the constraints 
of spacetime and internal symmetries. Starting from a 
classification of all operators allowed in the Lagrangian 
by symmetries, we propose two modifications: 

• Incorporation of explicit chiral symmetry breaking 
in the four-fermion interaction. 

• Addition of a two-derivative kinetic term, mimick¬ 
ing the lattice Wilson term. 

We then show that the modification of the four-fermion 
interaction is required for the model predictions to be 
consistent with lattice results, at least away from the 
continuum limit (at fixed lattice spacing). 

We would like to stress that it is not our purpose here 
to carry out a precise numerical fit of the effective model 
to available lattice data. We rather wish to gain qualita¬ 
tive insight with a reasonable number of free parameters, 
and thus to prepare the ground for a future more de¬ 
tailed quantitative study. The simple setting used here 
allows us to study separately, and make robust conclu¬ 
sions about, the various physical ingredients entering the 
problem, namely the chiral restoration at high density 
and low temperature as well as deconfinement at high 
temperature and moderate-to-high density. 

The outline of the paper is as follows. In section H 
we review the most important properties of 2cQCD (see 
Ref. [II] for a recent review). Section HI is devoted to a 
detailed discussion of the model. We classify all operators 
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up to dimension six, allowed by the symmetries, and give 
a brief overview of the qualitative changes brought about 
by the non-standard operators. Since we modify one of 
the interaction terms, it is mandatory to first analyze 
how it affects the physics in the vacuum. This is done in 
section IV, where we also fix the parameters in the quark 
(NJL) sector of the model. The next two sections consti¬ 
tute the core of the paper, where we investigate various 
lattice artifacts in accordance with our main program. In 
section V, we work at zero temperature, allowing us to 
separate the physics of flavor symmetry from deconfine¬ 
ment issues and the Polyakov loop. We analyze in turn 
the effects of a diquark source, the dependence on the 
quark mass, and the modified interaction term. In sec¬ 
tion VI, we add the gauge sector to the model and study 
the thermodynamics at nonzero temperature; the renor¬ 
malization of the Polyakov loop is a new element here. 
Subsequently, we extract the chemical potential depen¬ 
dence of the temperature scale in the Polyakov-loop po¬ 
tential, which essentially determines the position of the 
deconfinement crossover. In section VII we summarize 
and conclude. Some calculation details can be found in 
appendix A, where we analyze in detail the consequences 
of the two-derivative (Wilson) kinetic term for quarks. 

II. TWO-COLOR QCD 

In this paper, 2cQCD means a non-Abelian gauge the¬ 
ory with the SU(2) gauge group and Nf degenerate fla¬ 
vors of fundamental quarks. Two-color QCD differs in 
many respects from real-world, three-color QCD, most 
of them stemming from the fact that the fundamental 
representation of SU(2) is pseudoreal. For instance, a 
color singlet can only be made out of an even number 
of quarks, hence baryons in 2cQCD are bosons. Conse¬ 
quently, dilute nuclear matter is expected to be formed 
by a Bose-Einstein condensate (BEC) of bosonic baryons 
rather than by a Eermi sea of nucleons. 

Next comes the question of the low-energy hadronic 
spectrum, which is of importance for the thermodynam¬ 
ics at low temperatures. The lowest-lying states in the 
spectrum are determined by the spontaneously broken 
flavor symmetry in the (2c)QCD vacuum. Here the pseu¬ 
doreality of quarks implies that the usual SU(V/)l x 
SU(V;)r X U(1)b chiral symmetry of QCD in the limit 
of vanishing quark masses is embedded in an extended 
SU(2V/) flavor symmetry group [12]. The chiral con¬ 
densate in the vacuum breaks this to Sp(2V/), lead¬ 
ing to 2V| — Nf — 1 Goldstone bosons. These include 
Nj — 1 pseudoscalar mesons, and Nf{Nf — l)/2 diquark- 
antidiquark pairs. Nonzero (degenerate) quark masses 
make these modes massive. Eorming an irreducible mul¬ 
tiple! of the Sp(2V/) symmetry, these states are all de¬ 
generate with a common mass that we denote by 7Ti,r- 

The determinant of the Dirac operator of 2cQCD is 
necessarily real and for an even number of degenerate 
quark flavors it is also positive [13]. Consequently, the 


theory does not suffer from the sign problem and Monte 
Carlo simulations of dense matter are possible. In the 
following, we will focus exclusively on the simplest case 
of two quark flavors. In this case, the usual pion triplet 
is augmented by a single diquark-antidiquark pair carry¬ 
ing baryon number but no isospin. This determines the 
basic topology of the phase diagram of 2cQCD. Nonzero 
baryon chemical potential /tb breaks the flavor symme¬ 
try down to the usual SU(2 )l x SU(2)r x U(1)b chiral 
group; this is natural as the additional symmetry gener¬ 
ators following from the pseudoreality of the quark rep¬ 
resentation do not commute with the baryon number op¬ 
erator. For nonzero (degenerate) quark masses, only the 
SU(2)v X U(1)b subgroup is exact. 

When /tb > w,r, the diquarks are expected to undergo 
BEC, which breaks U(1)b spontaneously. Since U(1)b is 
an exact symmetry, the baryon superfluid phase is nec¬ 
essarily separated from the vacuum by a phase transi¬ 
tion. As the chemical potential is further increased, one 
eventually enters a Bardeen-Cooper-Schrieffer-like (BCS- 
like) regime where the thermodynamics is dominated by 
a Fermi sea of quarks that form weakly bound Cooper 
pairs. The order parameter for U(1)b symmetry breaking 
now is a composite diquark operator, which has the same 
quantum numbers as the order parameter in the BEC 
phase. The two regimes should therefore be smoothly 
connected [14]; in this context one speaks of a BEC-BCS 
crossover [15]. 

If we instead crank up the temperature at zero baryon 
chemical potential, we expect the physics to be more- 
or-less similar to that of three-color QCD. Around 
some pseudocritical temperature Tc, we expect a rapid 
crossover from hadronic to quark degrees of freedom, 
accompanied by a rise in the expectation value of the 
Polyakov loop. This is loosely referred to as deconfine¬ 
ment. In the same range of temperatures, the chiral 
condensate melts, and we enter the quark-gluon plasma 
phase. 

In lattice simulations of 2cQCD with Wilson fermions, 
the symmetry of the theory is affected in a twofold man¬ 
ner. Eirstly, Lorentz invariance is broken to a discrete 
symmetry group of the spacetime lattice, which may re¬ 
sult in a certain degree of anisotropy. At the same time, 
the Wilson term in the action acts as an additional source 
of explicit breaking of flavor symmetry which only disap¬ 
pears in the continuum limit. Since it breaks the flavor 
symmetry in exactly the same way as the quark masses 
do it may not he possible to distinguish the two sources 
of symmetry breaking in low-energy observables. 

In our model calculations, we will focus on two classes 
of observables: symmetry-breaking order parameters 
(condensates) and thermodynamic quantities such as 
pressure or baryon density. We now summarize the ex¬ 
pectations for these observables and then describe the 
relevant lattice results in order to set the stage for our 
analysis. 
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A. Chiral perturbation theory 


angle of rotation 9 is given by [12] 


The flavor symmetry of 2cQCD and its spontaneous 
breaking are most conveniently encoded in the low-energy 
effective theory for the pseudo-Goldstone modes, namely 
chiral perturbation theory (yPT) [16]. Following the no¬ 
tation of Refs. [12, 17], the pions and diquarks are ex¬ 
pressed in terms of a 2iVj x 2Nf antisymmetric unimodu- 
lar unitary matrix S, in terms of which the leading-order 
yPT Lagrangian in Minkowski space reads 

■^xPT = +i7Retr(JS). (1) 


Here is the pion decay constant and the covariant 
derivative includes the baryon number chemical po¬ 
tential -1- YB'^), where B 

is the baryon number operator. The matrix J in general 
contains sources that couple to scalar or pseudoscalar 
quark bilinears. Here we will need the scalar source toq 
( quark mass), and the diquark source j, in terms of which 
we have J = moY\ -|-jE^. In the case of two flavors, the 
matrix basis can be chosen such that 




Si 




in the 2x2 block form, where T 2 is a Pauli matrix. Fi¬ 
nally, the coupling H can be fixed by expanding the La¬ 
grangian (1) to second order in the fluctuations about 
the vacuum, or by using the Gell-Mann-Oakes-Renner 
relation, leading to H — 

As the low-energy spectrum contains excitations car¬ 
rying baryon number, the leading-order chiral La¬ 
grangian (1) can be used to study the phase diagram 
of 2cQGD at zero temperature and nonzero baryon den¬ 
sity [18]. The ground state at hnite density is most 
easily visualized by using a Lie algebra isomorphism 
to cast the coset space SU(4)/Sp(4) equivalently as 
SO(6)/SO(5) [17]. The unitary matrix E can thus be 
mapped onto a unit 6-vector n via the relation Y = ft -Y, 
where E^ is a set of suitably chosen basis matrices. De¬ 
pending on the values of the chemical potential(s) and 
the external sources, the ground state therefore moves 
on a unit sphere. In the absence of an isospin chemical 
potential fj-i and other sources except mo,j, the symme¬ 
try of the problem can be exploited to rotate the ground 
state into the ( 711 , 712 ) plane, corresponding to a chiral 
condensate and a diquark condensate with a fixed phase. 
The ground state can thus be parametrized by a single 
angle 9 such that tii = cos 9 and 712 = sm9. 

For the time being, we will assume that j = 0. The 
effects of the diquark source will be discussed in detail in 
section V A. The static part of the Lagrangian (1), whose 
maximum is to be found, then becomes 


sin^ 9 + 4/^771^ cos 9. (3) 

As expected, the chiral-symmetry-breaking state 0 = 0 is 
stable for /rp < Tn^r. As /tb further increases, the equilib¬ 
rium starts rotating into the diquark direction, and the 
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77? 

cos9=^. (4) 

This result is a priori expected to hold only within the 
range of validity of yPT, in particular only in the BEG 
regime where bosonic degrees of freedom dominate the 
physics of dense two-color matter. It should therefore be 
emphasized that, in fact, it remains at least qualitatively 
accurate even for much higher values of ^b- A numeri¬ 
cal solution of the NJL model shows that Eq. (4) holds 
also for chemical potentials rather deep in the BGS phase 
where a Fermi sea of quarks has been formed. The agree¬ 
ment between the two approaches, based on completely 
different degrees of freedom, persists at zero temperature 
up to /7 b ~ [19]. This can be attributed to the fact 

that in this range of fj-B , the physics at zero temperature 
is almost entirely driven by the condensates; the contri¬ 
bution of quark quasiparticles is negligible. 


B. The puzzle and its resolution 

The above discussion hints that the quantitative as¬ 
pects of the phase diagram of 2cQCD are determined 
by its symmetries even far beyond the region where one 
would expect it. The results of lattice simulations are in 
a stark contrast to this naive expectation. Most impor¬ 
tantly, the lattice data indicate a fast transition to a BGS 
regime just above /ts = the expected bosonic BEG 
phase, if present at all, is not resolved [20]. Moreover, in 
the BGS regime, 2cQCD behaves as a system of weakly 
interacting, almost massless quarks. This conclusion is 
supported by two independent pieces of evidence [3]: 

• The values of pressure and baryon number density, 
when normalized to their values for an ideal gas 
of massless quarks — the Stefan-Boltzmann (SB) 
limit — exhibit a plateau at /7 b ^ tti^ close to 
one. The precise height of this plateau is currently 
hard to determine, but its existence seems to be 
confirmed. 

• The expectation value of the color-singlet diquark 
operator scales with /tg in the same range of chem¬ 
ical potentials, reminiscent of the density of states 
at a sharp Fermi surface of massless relativistic 
fermions. 

An additional, related piece of evidence is provided by 
the fact that the critical temperature for diquark con¬ 
densation saturates at high /7 b at a value roughly given 
by « Tc/2 [4], whereas the picture of the order pa¬ 
rameter rotating on a unit sphere, sketched above, would 
naively suggest a vastly different value, « Tc. 

How can we reconcile these observations with the uni¬ 
versal model-independent predictions of xPT? The sim¬ 
ple and short answer is: we cannot. The SB scaling of 


4 


thermodynamic observables requires that the quarks are 
almost massless and their quasiparticle gap, proportional 
to the diquark condensate, is very small as well. It is ob¬ 
vious that this cannot be true simultaneously in yPT, 
according to which the sum of squares of the chiral and 
diquark condensates is constant. 

It is apparent that the key physical ingredient required 
is the rapid transition to a BCS-like gas of almost mass¬ 
less and gapless quarks at /ie ^ Once this is 

achieved, the other features — ^g-scaling of the expec¬ 
tation value of the diquark operator and the suppression 
of Til — should follow naturally. yPT as well as effective 
models based on the same symmetries lead to predictions 
that are at odds with this requirement. We therefore ex¬ 
pect explicit chiral symmetry breaking to play a crucial 
role. However, as we demonstrate in section VB, tun¬ 
ing the current quark mass even to unreasonably high 
values is not sufficient. In the remainder of this paper, 
we therefore take a rather radical approach to the prob¬ 
lem in which we abandon the chiral symmetry altogether. 
We introduce explicit symmetry breaking into the effec¬ 
tive four-quark interaction and show that this leads to 
the desired effect. The possible origin of this symmetry 
breaking is discussed in section VII. 


III. THE MODEL 


Our model is of PNJL type, which is based on quark 
degrees of freedom. Thus the only field variable is the 
quark spinor ijj. For the time being, we focus on the 
quark sector. The gauge sector will be discussed in de¬ 
tail later, when the effects of nonzero temperature are 
introduced. The form of the Lagrangian is constrained 
by spacetime and internal symmetries, which we assume 
to be as follows: Poincare invariance plus the discrete 
symmetries of charge conjugation, parity, and time re¬ 
versal, global SU(2) color symmetry, and Sp(4) flavor 
symmetry. We therefore abandon the full SU(4) flavor 
group. This group is explicitly broken down to the Sp(4) 
subgroup by (degenerate) quark masses and the Wilson 
term in the lattice action. The Sp(4) subgroup therefore 
constitutes the true flavor symmetry of both lattice and 
continuum 2cQCD with massive quarks. 

Finally, we make one more step which goes beyond 
the usual PNJL model building. We do not restrict our¬ 
selves to the simplest possible Lagrangian consisting of 
a quark kinetic term and a four-quark interaction. In¬ 
stead, we classify all terms in the Lagrangian consistent 
with the above symmetries up to dimension six [21]. This 
is in agreement with the effective field theory philosophy 
where all operators up to a given dimension allowed by 
symmetry should be included [22]. 


A. Classification of operators 

We start by classifying the operators according to their 
canonical dimension. The operators are written schemat¬ 
ically with their indices suppressed. 

• Order 3. The only parity-even Lorentz scalar that 
respects baryon number, SU(2) color and isospin 
invariance is the quark mass term tpip. The Sp(4) 
symmetry is automatically implied, but the full 
SU(4) group is explicitly broken. 

• Order 4. Schematically, the operator must be of 
the ipDip type to respect baryon number. Lorentz 
invariance and parity together with the SU(2) color 
and isospin invariance single out the usual quark 
kinetic term {plpij; . The full SU(4) symmetry is 
automatically implied. 

• Order 5. Here we have two possibilities respecting 

baryon number conservation and Lorentz invari¬ 
ance, namely and The 

latter is, however, irrelevant for our purposes since 
we only consider a purely temporal background 
gauge field, representing the baryon chemical po¬ 
tential and the Polyakov loop. Parity, SU(2) color 
and isospin invariance then single out the operator 

The Sp(4) symmetry is automatically 
implied, but the full SU(4) group is explicitly bro¬ 
ken. We will refer to this operator as the Wilson 
term since it closely resembles the corresponding 
operator in the lattice Wilson action. 

• Order 6. Here we have two schematic possibilities 
respecting baryon number, namely tpDDDip and 

. Lorentz invariance requires the former to 
contain an odd number of Dirac matrices, hence it 
represents a higher-order correction to the kinetic 
term which automatically preserves chiral symme¬ 
try. We will therefore drop it from further con¬ 
sideration since we do not expect that it leads to 
any qualitatively new effects. The non-derivative 
operator, on the other hand, is a standard NJL- 
type four-quark interaction. There are many oper¬ 
ators of this type that respect a given symmetry, 
as can be seen by performing a Fierz transforma¬ 
tion of the basic one-gluon-exchange type of oper¬ 
ator [7]. In the mean-field approximation, we want 
to build an invariant interaction out of fermion bi- 
linears that carry the quantum numbers of the low- 
energy degrees of freedom. The degrees of free¬ 
dom that must be present in the model are the 
pseudo-Goldstone bosons: the pion triplet (repre¬ 
sented by and the isospin-singlet diquark, 

which is represented by J5(J2T2'4> (where C stands 
for charge conjugation and CT 2 is the Pauli matrix 
in color space). Furthermore, the operator with 
the quantum numbers of a true scalar, must 
be added in order to account for the chiral conden- 
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sate. The two Sp(4)-invariant interactions that can 
be built out of the squares of these operators are 

and ('0i75r'0)^ + |V'‘^75cr2T2V'P- (5) 

The above considerations suggest the following generic 
NJL-type Lagrangian, 

^NJL = - mo)^ + + Gi'ip'ijjf 

- - ( 6 ) 

+ AG [(V’i75'rV')^ + , 

where the covariant derivative includes the baryon 
chemical potential, and later in section VI also the con¬ 
stant background gauge field representing the Polyakov 
loop. The dimensionless parameter A is in the remainder 
of this paper referred to as the chiral twist. The min¬ 
imal NJL model for three-color QCD amounts to skip¬ 
ping the last operator inside the brackets, and setting 
K = 0 and A = 1 [0]. The two-color version including 
the diquark channel, with k = 0 and A = 1, was first 
introduced in Ref. [23] and subsequently used by several 
other groups [24, 25]. In the following subsection, we will 
discuss the consequences of the modifications introduced 
here. 


B. Effects of chiral symmetry breaking 

The classification of invariant operators has automati¬ 
cally guided us to introduce two new couplings in the La¬ 
grangian. Let us start with the Wilson term proportional 
to K. As already mentioned, it mimics the discretization 
artifacts introduced by the lattice Wilson action. The 
presence of a two-derivative bilinear operator in the La¬ 
grangian leads to doubling of the fermion degrees of free¬ 
dom. For small n, we expect the new fermion species 
to be heavy with mass scaling as 1/k. Since the Wilson 
term breaks chiral symmetry explicitly, it will enhance 
the chiral condensate in the vacuum. Regarding thermo¬ 
dynamics, we expect the effects of the Wilson term to 
be most pronounced at high temperatures since the new 
fermion is heavy and thus difficult to excite thermally. 
In addition, the presence of a heavy fermion may lead to 
undesirable artifacts at high baryon density, namely the 
appearance of a second Fermi surface at high p,B and low 
temperature. Some of the quantitative consequences of 
the Wilson term will be worked out in appendix A. 

It is the chiral twist that will play a pivotal role in 
our analysis. As we abandon chiral symmetry here, the 
scalar and pseudoscalar channels are no longer forced to 
appear with the same strength. However, the diquark 
and pion operators do enter through a fixed combination, 
as required by the exact Sp(4) symmetry of two-flavor 
2cQCD. This guarantees exact mass degeneracy of pions 
and diquarks, and therefore that the critical chemical 
potential for diquark condensation at zero temperature 
equals the pion mass. 


A detailed analysis of the consequences of our 
model (6) with the modified interaction constitutes the 
bulk of the remainder of the paper. Yet some qualita¬ 
tive observations can be made already now. It is clear 
that the existence of a bound state in the pion channel 
requires A > 0. Since the strength of attraction in the 
scalar channel remains fixed to G, the constituent quark 
mass in the vacuum will be unaffected by A. As a conse¬ 
quence, the pion will become less strongly bound and its 
mass will increase with decreasing A. The location of the 
BEC-BCS crossover at zero temperature will therefore 
shift towards lower values of pB /Wtt and it is conceivable 
that the BEG region can be eliminated altogether. In 
the following, we will therefore consider only the range 
0 < A < 1. 

Since the coupling in the diquark channel also de¬ 
creases with A, we expect the diquark condensate to be 
suppressed. This will sharpen the quark Eermi surface 
and move the system closer to the SB limit. An addi¬ 
tional effect which follows from our analysis below, is that 
the chiral condensate will also be suppressed as compared 
to the naive expectation based on Eq. (4). Suppressing 
both the diquark and the chiral condensate, the chiral 
twist is therefore exactly what we need in order to ex¬ 
plain the rapid crossover to the gas of weakly coupled 
almost massless quarks. 

Erom the strict symmetry point of view, the Wilson 
term and the chiral twist should be treated on the same 
footing as they break symmetry in the same way. How¬ 
ever, a detailed numerical investigation reveals that the 
Wilson term has a rather marginal effect on the observ¬ 
ables considered in this paper, and certainly cannot re¬ 
solve the puzzle presented in section HB. We therefore 
take the liberty and discard the Wilson term altogether 
from most of our analysis. This step is driven by simplic¬ 
ity: as explained in detail in appendix A, including the 
Wilson term requires a substantial modification of the 
framework including a change of regularization scheme. 
At the level of accuracy, considered in this paper, the chi¬ 
ral twist turns out to be the only necessary and sufficient 
new ingredient in the model. 


IV. VACUUM PHYSICS 

In the remainder of the paper, except for appendix A, 
we set K = 0. The composite bosonic degrees of 
freedom cr, tt. A, and A* are introduced by the 
Hubbard-Stratonovich transformation, adding to the La¬ 
grangian (6) the term 

AifNjL = - + 2AGVJi75T'(/')^ 

- - 2AG^^i75cr2T2V’^l^. (7) 

The Lagrangian is now bilinear in the quark fields and we 
can integrate them out exactly. This leads to the effective 
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FIG. 1. Pion mass (left panel) and decay constant (right panel) as a function of A. All other parameters are fixed according to 
Eq. (19). 


bosonic action 


where the quark dispersion relation in the vacuum is 




off 

NJL 



4G 4AG J 


-iTrlog^, ( 8 ) 


Cfe = y/k^ + M^. (12) 


where 


, _ f 10 — M — iqsTT • T 

-A*75 


A 75 

10 — M + 1757 ? • T 


(9) 


is the Dirac operator acting on the Nambu space spanned 
by red quarks and green antiquarks. We have introduced 
the shorthand notation M = toq + cr for the constituent 
quark mass, and denoted by “Tr” the functional trace. 

In the following, we work in the mean-field approxi¬ 
mation. This means that after taking functional deriva¬ 
tives of the action ( 8 ) as appropriate, all bosonic fields 
are set equal to a constant. The integrals that appear 
are ultraviolet divergent and we need to regulate them. 
With the exception of appendix A, we will use a sim¬ 
ple sharp three-momentum cutoff A. Depending on the 
cutoff, writing the expressions in a manifestly Lorentz- 
covariant form may involve some manipulations that are 
not justified from a strictly mathematical point of view. 
This is, however, a well-known issue [ 6 ]. 

The vacuum of 2cQCD is characterized by broken chi¬ 
ral symmetry. Differentiating Eq. ( 8 ) with respect to mg, 
we find that the chiral condensate of a single quark flavor 
is related to the cr condensate by 

The latter is found from the stationary point condition 
= 0. This yields [26] 


The spectrum of pseudo-Goldstone bosons can be deter¬ 
mined from the polarization function (inverse propaga¬ 
tor), which is obtained by taking a second functional 
derivative of the action. We already know that in the 
vacuum pions and diquarks are degenerate, so we just 
state the result for the pion polarization function, sim¬ 
plified by using Eq. (11), 


(1 ■ (13) 


where 



10k 1 

(27r)4 [{k -f p)2 - M2](fc2 - M2) ■ 


(14) 


In the chiral limit, M = a and we can immediately see 
the effect of explicit chiral symmetry breaking by the 
chiral twist: the pion is exactly massless only for A = 1. 
In general the pion mass squared is given by the zero of 
the inverse propagator, x(m2) = 0. Using the the gap 
equation ( 11 ) once more, we can rewrite this as 


1 = 8AGA^c 


d^fc / 1 1 \ 

(27r)3 2 efc -I- 2tk - J ' 


(15) 


Using the Lehmann spectral representation of the pion 
propagator, we can determine the coupling g^gg of the 
one-pion state to the pseudoscalar quark bilinear 


cr = 16iGWM 


0k 1 

0100 00^00 


= SGNcM 


0k 1 

( 11 ) 


1 

Qirqq 


x'{p^) 


p^—m^ 


167V, 


d^fc Cfc 

( 16 ) 
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From this result, one can obtain the coupling of the one- 
pion state to the axial vector current, that is, the pion 
decay constant /,r- The resulting expression is 


u 


dirgg f ^ \ 

2Gml \\ )■ 


(17) 


A. Parameter fixing 

The NJL model (6) with k = 0 and A = 1 is defined 
by the three parameters mo, G, and A, which should 
be determined by a fit to three independent observables. 
It is customary to use the chiral condensate, pion de¬ 
cay constant and pion mass for that purpose. We follow 
Ref. [24] and determine the values of these input quan¬ 
tities in 2cQCD from their physical, three-color counter¬ 
parts using a naive scaling with the number of colors 

(uu) = -(218 MeV)^ 

= 75.4 MeV, (physical input) (18) 
rriTr = 140 MeV. 

Equations (10), (15), and (17) then give us the following 
values 

G = 7.23 GeV■^ 

A = 657 MeV, (fitted parameters) (19) 
mo = 5.4 MeV, 

which we will use throughout the rest of the paper un¬ 
less explicitly stated otherwise. The coupling A will be 
treated as a tunable parameter. 


B. Role of the chiral twist 

As explained in section IIIB, we expect that the pion 
mass increases with decreasing A. It is clear from the left 
panel of figure 1 that the effect is actually rather large. 
At A « 0.95, the pion mass is increased by a factor of 
two, and at around A « 0.6, it reaches the (unphysical) 
threshold for decay into a quark-antiquark pair. (The 
a condensate in the vacuum is to great precision equal 
to 300 MeV for our choice of parameters.) At the same 
time, the pion decay constant starts to drop rapidly (right 
panel of figure 1), which indicates that the pion ceases 
to behave as a Goldstone boson. The interval [0.6,1] 
therefore defines the range of reasonable values for A to 
which we will restrict ourselves in the following. 

To gain a more analytic insight into the A-dependence 
of physical observables, we take the derivative of Eq. (15) 
with respect to A. Using Eq. (16) allows us to rewrite the 
result in the form of an exact differential equation. 


dm^ 


^Trqg 

'2A2G' 


( 20 ) 


In the chiral limit, the pion mass will scale asymptotically 
as m^ oc Vl — A for A —>■ 1. Using in addition Eq. (17), 


we readily obtain another differential equation, this time 
for the pion decay constant. 


= 
J TT 


1 


{M - Xa)^ 


d 1 
dA m2 


( 21 ) 


This expression together with the asymptotic scaling 
mTT oc Vl ~ A explains the very weak A-dependence of 
f-jr for A close to one, see the right panel of figure 1. The 
fact that is almost constant in the range A S [0.8,1] 
also justifies a posteriori our treatment of A as a tunable 
parameter: since the quark mass, and hence the pion 
mass, is a free parameter on the lattice anyway, we have 
some freedom in tuning both mo and A without affecting 
the physical observables (uu) and and do not have to 
refit our parameters anew for each value of A. 


V. ZERO TEMPERATURE: CHIRAL 
RESTORATION 


In this section, we will investigate the effects of tun¬ 
ing various parameters on different physical quantities at 
zero temperature. We focus on the chiral and diquark 
condensates, the pressure, and the baryon number den¬ 
sity. We can therefore drop the pion field tt in Eq. (8). 
Evaluating the action for constant a and A and going to 
Euclidean space gives the thermodynamic potential 


,T= 

'NJL 


4G 


4AG 


± 


d^k 

(27r)' 


E 


k ’ 


( 22 ) 


where the quark quasiparticle dispersion relations are de¬ 
fined by [27] 


= v/(^^)" + |AP, = Cfc ± E (23) 

The values of the condensates for a given ^ are found 
by direct numerical minimization of the thermodynamic 
potential. The pressure is equal to —convention¬ 
ally shifted by a constant so that the pressure is zero in 
the vacuum. The baryon number density is obtained by 
taking the derivative of the pressure with respect to /jbj 
giving 


n 


T=o 

B 


= N, 


r d^k 



1 (27r)3 




(24) 


where we have separated the particle and antiparticle 
contributions. 


A. Role of the diquark source 

The first lattice simulations of dense 2cQCD with Wil¬ 
son quarks were performed with a fixed external source 
j for the diquark operator [28]; only recently has the ex¬ 
trapolation to vanishing source been studied. However, 











1.0 



j 

FIG. 2. Dependence of the rescaled diquark condensate <5 on 
the dimensionless source j for several different values of the 
parameter x = (given in bold). 


these attempts were based on fitting three data points by 
a simple analytical ansatz for the j-dependence. Ref. [3] 
resorted to a linear extrapolation, but as pointed out in 
Ref. [4], none of the three ansatze used therein (linear, 
power-law, and power-law with an offset) led to satisfac¬ 
tory results. 

In principle, adding a diquark source to the NJL model 
is straightforward and is analogous to introducing the 
quark mass toq as a source for the '0'0 operator. All 
we need to do is to make the shift A —^ ^ + j in the 
fermion part of the thermodynamic potential (22). How¬ 
ever, since we want to gain analytic insight into the scal¬ 
ing of the A condensate in the limit j —> 0, numerical 
solution of the NJL model is not satisfactory. Instead 
we once again employ yPT. We have confidence in this 
since the NJL model and yPT give numerically very sim¬ 
ilar results in a large range of chemical potentials [19]. 

We showed below Eq. (1) how the diquark source enters 
the yPT Lagrangian. Following Ref. [12], we relate the 
diquark source to the quark mass by introducing a new 
angle 4>, and express the chemical potential in terms of a 
dimensionless parameter x, via 

j = motan(j) = moj, x=^^. (25) 

The static part of the Lagrangian (3) with the added 
diquark source term can then be rewritten as a dimen¬ 
sionless potential, 


Vie) ^ - 


rc^stat 

‘=^XPT 

4/2to2 


„ „ = —-a: sin 0 — COS0 — Jsin0 

Zry-y-)Z Q 


^ 2 ■ 2n COS(0 - 4>) 

--x^sin'^e -^ 

2 cos(p 


(26) 


Finally introducing the shorthand notation 6 = sin 9 for 
the normalized diquark condensate, the stationarity con¬ 
dition 6V{9)/59 = 0 takes the simple form, 

j = ^ — x^S. (27) 


Note that for j = 0, we recover the nontrivial solution 
for the chiral condensate (4), now expressed as cos 9 = 
\/l — = 1/a;^. It is easy to see that for any j > 0, 

Eq. (27) admits a unique positive solution 6{j). 

The presence of the condensate in the ground state is 
reflected by a specific asymptotic scaling of this solution 
in the limit j —> 0. The asymptotic expansion of J(J) can 
be found by an iterative solution of Eq. (27), leading to 


S{j) 

S{j) 

SO) 


J 


1 — a:^ 2(1 — x'^)^ 


(2jr 


-b 


24 

j 


+ 0(?), 

Oi~f), 


x* x'^ix'^ — 1) 


(x < 1), 

{x = 1), 


(28) 


3x0^ 

2(a;4 _ 1)5/2 


Oi'f), 


{x > 1). 


The a; = 1 part is most easily obtained by writing 
Eq. (27) as j = 4tan^ |/(1 — tan"^ |), first solving for 
9ij), and then converting this into a series for J(J). In 
contrast, the solutions in the regions a: < 1 and a: > 1 
are simple Taylor expansions. Moreover, for a: < 1 the 
expansion only contains odd powers of j, reflecting the 
unbroken discrete symmetry of Eq. (27), under which 
both j and 6 change sign. We stress that in both regions 
a: < 1 and a: > 1, the condensate is a linear function of 
the source in the limit j —>■ 0, but the series convergence 
becomes slower and slower as the phase transition is ap¬ 
proached. A rough upper bound on the range of values 
of j in which a linear extrapolation makes sense, can be 
obtained by comparing the first two j-dependent terms 
of the expansion, giving j < |a^ — for x close to 
one in both regions. The convergence of the expansion 
can also be judged from the exact numerical solution of 
Eq. (27), shown in figure 2. 


B. Role of the quark mass 

We argued in section IIB that reproducing the fast 
transition to a gas of weakly interacting almost massless 
quarks at /te > m^n-, as seen on the lattice, requires strong 
explicit breaking of chiral symmetry. It might be tempt¬ 
ing to the think that this effect is due to a large current 
quark mass. In this section, we therefore set A = 1 and 
investigate the dependence on mg at zero temperature 
using the mean-field approximation (22). 

Since the transition to the baryon superfluid phase is 
expected to occur at /tb = m-Tr, it makes sense to trade 
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FIG. 3. Dependence of the rescaled chiral and diqnark condensate on a: = for several different values of the current 

quark mass mo- Here ctq is the chiral condensate in the vacuum for given mg and the other parameters as in Eq. (19). The 
values of mg in MeV are indicated in bold (the same color coding is used in both panels). The dashed line shows the prediction 
of xPT. 
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FIG. 4. Dependence of the rescaled chiral (left panel) and diquark (right panel) condensate on x = for several values 

of A (shown in bold). Here ag « 300 MeV is the chiral condensate in the vacuum for our parameter set (19). The dashed line 
indicates the prediction of xPT. 


the chemical potential for the dimensionless parameter 
X = for the sake of comparison. This is done 

in figure 3, where the solutions to the gap equations for 
the chiral and diquark condensates are shown for sev¬ 
eral values of mo- It is obvious, after proper rescaling, 
that the prediction of xPT (4), indicated by a dashed 
line in the figure, works quite well even for unreasonably 
heavy quarks. (One would probably not expect an ap¬ 
proach based on spontaneously broken symmetry to still 
work even when the mass of the pseudo-Goldstone bo¬ 
son reaches the scale of the ultraviolet cutoff.) It is also 
clear that just tuning the quark mass not is sufficient for 


our purposes, even for the largest values of mo: the size 
of the quasiquark gap A remains largely unaffected [29]. 
Moreover the slight reduction in the chiral condensate is 
canceled by the increased current quark mass so that the 
constituent quark mass is not reduced at all. 

The results shown in figure 3 are also in a good 
agreement with older lattice simulations using staggered 
quarks [13, 30]. There, the prediction (4) of yPT was ver¬ 
ified numerically for current quark masses varying by an 
order of magnitude. A tiny reduction of the chiral con¬ 
densate as compared to Eq. (4), observed therein, could 
— exactly as in our case — be ascribed to the relatively 
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FIG. 5. Dependence of baryon number density (left panel) and pressure (right panel) on chemical potential for several values 
of A (shown in bold, the same color coding as in figure 4 is used). Both quantities are normalized to the SB values for a gas of 
free massless quarks. 


large current quark mass. Since the staggered implemen¬ 
tation of lattice quarks preserves chiral symmetry, we 
conclude that heavy quarks alone are not sufficient to 
explain the thermodynamic behavior in the baryon su¬ 
perfluid phase observed in Ref. [3]. An additional source 
of chiral symmetry breaking is needed. On the lattice, 
this is provided by the Wilson term in the action. In our 
NJL model, the chiral twist serves this purpose, as we 
will now demonstrate. 


C. Role of the chiral twist 

In this section, we fix mo to its value given in Eq. (19), 
and instead vary the chiral twist. Figure 4 shows the de¬ 
pendence of the chiral and diquark condensates at zero 
temperature on A. Since the pion mass is very sensi¬ 
tive to A, we again plot these quantities against the ratio 
X = The numerical results fully confirm our 

expectations outlined in section IIIB, namely that the 
diquark condensate (at fixed x) becomes rather strongly 
suppressed as A decreases. At the same time, the chiral 
condensate is strongly suppressed as well. We note that 
this is in contrast to yPT where the sum of the two con¬ 
densates squared is constant. The suppression of both 
condensates simultaneously is exactly what we want. 

How small should A be to reach a quantitative agree¬ 
ment with the lattice results of Refs. [3, 4]? The sup¬ 
pression of the diquark condensate leads to a suppres¬ 
sion of the critical temperature Td- In this paper, we are 
not going to compute the critical temperature, but we 
can make at least an estimate [31]. Since in the weak- 
coupling BCS theory as well as in its extension including 
the most important corrections due to fluctuations, the 


critical temperature is proportional to the pairing gap at 
zero temperature, and since we want Td ~ T'c/2, we need 
to reduce the diquark condensate roughly by a factor of 
two. This suggests that A should fall somewhere in the 
range [0.6, 0.7]. 

Let us next consider the thermodynamic observables 
baryon number density and pressure. (Since we are at 
zero temperature, the energy density is linearly depen¬ 
dent on these two via the Gibbs-Duhem relation.) Both 
quantities are plotted in figure 5. The peak just after the 
onset of diquark condensation at /tb = well formed 
at A = 1, is a hallmark of BEG. As A decreases, it is 
gradually washed away, and again we need A G [0.6,0.7] 
to be able to conclude that the BEG-like behavior has 
given way to BGS-like scaling. One should, however, be 
aware of the fact that the absolute height of the peak is 
only indicative here, as it is very sensitive to the value of 
the pion decay constant and thus to the precise values of 
the parameters used [19]. 

Before closing this section, we would like to point out 
a subtle detail concerning the dependence of our results 
on the chiral twist. As may be seen already in figure 4, 
the diquark condensation phase transition becomes first 
order for sufficiently small A. In order to highlight this 
feature, we zoom in on figure 4 in the immediate vicin¬ 
ity of the transition. This is shown in figure 6. The 
discontinuity of the condensates becomes rather strong, 
yet the location of the transition remains very close to 
= 1. For the values of A considered here, it re¬ 
mains in the range of ^b/™ 7 i- S [0.98,1]. We have con¬ 
firmed numerically that around the transition point, the 
thermodynamic potential has two competing local min¬ 
ima, and used this to locate the transition precisely. The 
first-order transition only appears for sufficiently low val¬ 
ues of A, the critical value being approximately Ac ~ 0.88. 








11 


o 

b 




FIG. 6. Detail of figure 4 close to the phase transition. The same color coding is used. 



FIG. 7. Dependence of the rescaled chiral and diquark condensate on the diquark source j for A = 0.8. The values of j in MeV 
are indicated in bold. 


While the appearance of a first-order transition was 
somewhat unexpected to us, it is not in contradiction 
with symmetry or any other physical principle. The 
change of the order of the transition may well be an arti¬ 
fact of our model, and we therefore do not analyze it any 
further. It is nevertheless interesting to note that intro¬ 
ducing the chiral twist provides a mechanism for a direct 
transition from the vacuum to the BCS regime of weakly 
coupled quarks without the necessity for an intermediate 
BEG phase. In any case, we remark that the first-order 
transition is unlikely to be visible in current lattice sim¬ 
ulations, since it appears very close to its expected po¬ 
sition, and since even a small external diquark source j 
is likely to turn it into a crossover. This is clear from 
figure 7, showing details of the transition for A = 0.8: an 
external source as small as 1 MeV is sufficient to smooth 


out the transition into a crossover. 


VI. NONZERO TEMPERATURE: 

DECONFINEMENT 

At nonzero temperature, the nature of thermal excita¬ 
tions in the system becomes important, and the simple 
NJL model (6) does not capture the equilibrium ther¬ 
modynamics of (2c)QCD correctly: at low temperatures, 
the relevant degrees of freedom of 2cQCD are the pions 
and the diquarks, whereas the NJL model is based on 
the quark degrees of freedom. We follow the by now 
standard procedure and take into account the confining 
property of strong interactions by coupling the model to 
the Polyakov loop. We have avoided doing so until now 
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for two reasons: (i) there is considerable freedom in the 
choice of the gauge sector of the model [32], and (ii) the 
resulting framework is no longer a Lagrangian field the¬ 
ory in the usual sense but merely a statistical model, 
since there are no dynamical gauge degrees of freedom. 

At nonzero temperature, we have to add the effects of 
thermal quark excitations as well as the gauge sector to 
Eq. (22). The full thermodynamic potential of the PNJL 
model 2cQCD then becomes [24] 


^PNJL — ^gauge(^) 

d^fe r 


a 


4AG 


-2N, 


(29) 


E 


(27r)3 L 


-b riog(l -b 2$e-^^'= -b ) 


where $ is the expectation value of the Polyakov loop 
and ngauge(d>) is the yet unspecified contribution of the 
gauge sector. 


A. Renormalization of the Polyakov loop 


Before we proceed with the numerical solution of the 
model, we have to address a conceptual issue related to 
the Polyakov loop. On the lattice, the concept of the 
Polyakov loop is rather subtle as in the naive contin¬ 
uum limit, its expectation value vanishes even in the 
deconfined phase; the Polyakov loop requires renormal¬ 
ization [33]. A simple way to think of it is as follows. 
The Polyakov loop expectation value can be related to 
the free energy of a static heavy quark immersed in the 
colored medium, Fq, via $ ~ An additive renor¬ 

malization of the free energy then gives a multiplicative 
renormalization of the Polyakov loop via 

(30) 


where <i)o is its “bare” value. On the lattice, another in¬ 
terpretation of this formula is available. Using the rela¬ 
tion /3 = Nr as, where Nr is the number of lattice points 
in the temporal direction and ag is the lattice spacing, 
we have . The lattice 

Polyakov loop is a time-ordered product of Nr link vari¬ 
ables, winding around the temporal direction. Renormal¬ 
ization of each link variable by the constant factor 
therefore gives rise to a temperature-dependent renor¬ 
malization of the Polyakov loop. 

In practice, the Z,j, factor is found by imposing a cer¬ 
tain renormalization condition. Following Refs. [3, 4], we 
define this condition by prescribing a value for the 
renormalized Polyakov loop at a given reference temper¬ 
ature T and p-e = 0. This leads to the relation 


= ^o{T,hb) 




T/T 


L^o(T,0)J 


(31) 


The factor is kept constant throughout the compu¬ 
tation: it is fixed at /xb = 0 and the same value is used 
regardless of the baryon chemical potential. 



T [MeV] 

FIG. 8. Expectation value of the rescaled chiral condensate 
(black line) and the Polyakov loop (red line) at /xb = 0 as 
a fnnction of temperature. The data points are taken from 
Ref. [4]; the highest data point is the reference point which 
defines the renormalization condition. Note that in order to 
be able to display both curves on the same scale, the chiral 
condensate is rescaled to equal 1/2 in the vacuum. 


How should we compare the prediction of our 
model (29) to the renormalized Polyakov loop measured 
on the lattice? Note that depending on the reference 
value the renormalized Polyakov loop can in princi¬ 
ple take on any positive value, whereas the quantity <i> 
in Eq. (29) is usually assumed (and for some choices of 
^^gauge enforced) to lie in the range [0,1]. Here we take 
the point of view that in our model, the finite-valued 4> 
is already renormalized, but in some a priori unknown 
scheme. In order to be able to compare our model to 
lattice results, we therefore first have to minimize the 
thermodynamic potential (29) with respect to 4>, and 
then perform an additional finite renormalization using 
Eq. (31) with the same renormalization condition as used 
on the lattice. 


B. Zero chemical potential 

We next have to make a choice for the Polyakov loop 
potential Hgauge(‘f’)- Some of the potentials used in the 
literature include a large number of free parameters, al¬ 
lowing a precise numerical fit to lattice data [9]. However, 
since we aim at a qualitative understanding rather than 
numerical fitting, we prefer to have a model with as few 
free parameters as possible. We therefore employ the po¬ 
tential already used in the context of 2cQCD in Ref. [24] , 
motivated by the lattice strong-coupling expansion [8], 

f^gauge($) = -6T[24$2e-/3« + log(l - $2)], (32) 

The parameter a is proportional to the deconfinement 
temperature Tg in the pure-gauge theory via a = 
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FIG. 9. Expectation value of the renormalized Polyakov loop for A = 0.7 and several values of /is. The red lines represent the 
best fit of the Tg parameter. The orange bands were obtained by varying Tg in order to estimate the error of the fit. The data 
points are taken from Ref. [4]. 


Tg log 24. The parameter b can be related to the string 
tension. While the parameters of the quark sector of the 
model are fixed by Eq. (19) — the chiral twist does not 
affect the mean-field thermodynamic potential at /ip = 0 
— the parameters a, b will be determined by a fit to lat¬ 
tice data for the expectation value of the Polyakov loop. 

To this end, we use the conversion between the lattice 
and physical units provided by Ref. [3], defined by the 
reference temperature of f = 281 MeV for Nt = 4. We 
then use the renormalization scheme B of Ref. [4], in 
which = 1/2. Our best fit to the lattice data is shown 
in figure 8 and corresponds to the values 

b = (278 MeV)^ Tg = 247 MeV. (33) 

In the same figure, we also display the temperature de¬ 
pendence of the chiral condensate (rescaled for conve¬ 
nience to equal 1/2 in the vacuum). This demonstrates 
that with the present values of the gauge sector parame¬ 
ters, differing somewhat from those of Ref. [24], the chiral 


and deconfinement crossovers at zero chemical potential 
appear in the same range of temperatures. Note that the 
location of the crossover may depend somewhat on the 
choice of the renormalization condition. We prefer not to 
give a single value for a pseudocritical temperature due 
to the ambiguity of this concept. 


C. Chemical potential dependence 

With increasing chemical potential, one generally ex¬ 
pects that the deconfinement crossover moves towards 
lower temperatures due to the back-reaction of the dense 
medium to the gauge sector. However, it is notoriously 
difficult to take this back-reaction into account in the 
PNJL model. A common way around this is to include 
an explicit /re-dependence in the potential flgaugej usu¬ 
ally using an analytic ansatz motivated by perturbation 
theory [34, 35]. One should note that introducing such a 
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chemical potential dependence into Ugauge leads to an 
unphysical artifact: it gives an extra contribution to 
the baryon number density that does not arise from the 
quark degrees of freedom of the model. However, since 
f^gauge = 0 at T = 0, this artifact does not violate the 
“Silver Blaze” property stating that at zero temperature, 
physics is independent of /tb below the onset of diquark 
BEC [36]. 

Here we would like to make the case that the data for 
the expectation value of the Polyakov loop at nonzero 
baryon density of Refs. [3, 4], can be exploited to obtain 
a direct information about the Polyakov loop potential, 
without resorting to any ansatz. To make it as simple 
as possible, we assume that b is constant and fixed by 
Eq. (33), while the parameter a contains a /ie-dependent 
temperature scale, 

a(/XB) = Tg(^B)log24. (34) 

This is in accord with the ideas put forward in Ref. [34], 
based on the perturbative running with the physical scale 
set by the chemical potential. Figure 9 shows a compar¬ 
ison of our model predictions with lattice data for dif¬ 
ferent values of fiB taken from Ref. [4]. In each plot, we 
have adjusted the value of Tg{^B) to represent the lattice 
data most faithfully; the quality of the fit can be assessed 
by varying Tg, indicated by the bands in figure 9. Along 
with adjusting Tg{^B), the /iB-independent value A = 0.7 
was chosen to achieve the best overall fit. Note that this 
value agrees well with the estimate A S [0.6,0.7] obtained 
in the previous section, based on different observables for 
a different thermodynamic regime, namely high density 
and zero temperature. 

The extracted values of Tg{^B) are shown as the blue 
data points in figure 10, with error bars defined by the 
bands in figure 9. This plot provides a direct informa¬ 
tion about the back-reaction of the dense medium to the 
gauge sector, without the need to employ a particular 
analytical ansatz for the function Tg(fiB)- If desired, the 
data points in figure 10 can of course be fitted with a 
suitably chosen function of /iB. 

The effect of tuning the chiral twist on the Polyakov 
loop crossover can be appreciated with the help of fig¬ 
ure 11. For A = 1, the crossover tends to be too steep. 
The reason is that we are in the baryon superfiuid phase 
and the quark gap A is too large to allow thermal ex¬ 
citation of quarks. Hence the behavior of the Polyakov 
loop to a large extent is the same as in the pure gauge 
theory. Once A is lowered, the quark gap is reduced and 
the light quark excitations smear out the crossover. The 
fact that the curves for different A converge at high tem¬ 
peratures is easy to understand: at these temperatures, 
the system has already undergone the phase transition 
to the normal phase (visible in some of the curves as a 
small cusp) where the thermodynamic potential is inde¬ 
pendent of A altogether. The residual deviation between 
different curves is a consequence of rescaling /iB by the 
A-dependent pion mass for the sake of the plot. 

Before closing the section, two remarks are in order. 



/is/m^ 


FIG. 10. Temperature scale Tg of the Polyakov loop poten¬ 
tial (32) as a function of fislrnn- The blue data points are 
extracted directly from figure 9; the error bars are defined 
by the bands therein. The red data points correspond to 
mo = 48 MeV, with the other parameters of the NJL model 
still given by Eq. (19); this reproduces the pion mass used 
in Ref. [4]. Both data sets correspond to the same values of 
fJ^s/rnTr, but the red points were slightly displaced for the sake 
of convenience. 



FIG. 11. Expectation value of the renormalized Polyakov loop 
for = 1.24 and several different values of A. The same 

color coding as in figure 4 is used: decreasing A by steps of 
0.1 moves from red (A = 1) to cyan (A = 0.6). 


Firstly, we have not made an attempt to fit all the pa¬ 
rameters of our model precisely to the lattice data. As a 
consequence, our pion mass (about 550 MeV for A = 0.7) 
differs by about 25% from its lattice value (720 MeV in 
physical units). The latter can be reproduced in our 
model by setting mg = 48 MeV while keeping the other 
parameters in Eq. (19) unchanged. Following the same 
procedure, that is fitting the parameters a, 6 at /jb = 0 to 
the lattice data and then re-adjusting Tg for each partic- 
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ular value of fiB, we arrive at the red data points shown 
in figure 10. The proximity of the two data sets demon¬ 
strates that our results are rather robust and to a large 
extent independent of the precise model setup. Note that 
in this latter calculation, including a nonzero diquark 
source is essential to achieve a good agreement with the 
lattice data. The effect of the diquark source on the ex¬ 
pectation value of the Polyakov loop is much weaker for 
the light quark parameter set (19), and all results shown 
in figures 8 and 9 are understood at j = 0. 

Secondly, as can be seen in figure 9, our model cannot 
reproduce very well the Polyakov loop data at both low 
and high temperatures. Reaching agreement here may 
require going beyond the simple model of Eq. (32) and 
introducing a more phenomenological gauge sector po¬ 
tential with a larger number of free parameters, or using 
input from other methods [37]. Nevertheless, for lower 
temperatures and lower /tb, our simple model can repro¬ 
duce the lattice data quite well, which makes our conclu¬ 
sions based on symmetry and its explicit breaking rather 
robust. 


VII. CONCLUSIONS AND OUTLOOK 

In the present paper, we have studied in detail the 
PNJL model with focus on comparing it to lattice data 
for 2cQCD. We have discussed extensively the effects of 
strong explicit chiral symmetry breaking due to finite 
quark masses and the Wilson term in the lattice ac¬ 
tion. We argued that taking these effects into account 
requires a generalization of the standard PNJL model. 
This is accomplished by introducing the chiral twist that 
incorporates explicit breaking of chiral symmetry in the 
four-quark interaction term. Based on diverse pieces 
of evidence, we have argued for its value in the range 
A € [0.6, 0.7], which can simultaneously account for the 
following nontrivial effects: 

• Rapid crossover from the vacuum to a BCS-like gas 
of almost massless and gapless quarks at /tb ^ nT-n-- 

• Strong suppression of the critical temperature 
for baryon number breaking as compared to the 
(pseudo)critical temperature Tc for chiral restora¬ 
tion at p-B = 0. 

• Strong broadening of the Polyakov loop crossover 
within the baryon superfluid phase. 

We argued in section VB that these effects cannot be 
explained by a large current quark mass. However, 
this statement should be understood within the (P)NJL 
model framework. We cannot discriminate with confi¬ 
dence the effects of explicit chiral symmetry breaking 
by lattice discretization and by the current quark mass 
on the level of the 2cQCD Lagrangian. Which of these 
sources of symmetry breaking is more important can be 
decided with future lattice simulations by varying the 
quark mass and lattice spacing. 


In addition, we were able to extract the /iB-dependent 
temperature scale Tg{^B) from the lattice data for the 
expectation value of the Polyakov loop. This is a model- 
independent result in the sense that it does not rely on a 
particular analytical ansatz for the /^B-dependence of the 
Polyakov-loop potential. However, one should not inter¬ 
pret Tg as a deconfinement temperature: while the actual 
transition from hadronic to quark degrees of freedom is 
a smooth crossover which does not have a well-defined 
critical temperature, the quantity Tg is well defined, but 
depends on the choice of the Polyakov loop potential. 

Our simple mean-held model of course has its limita¬ 
tions, most importantly in that it ignores huctuations of 
the composite bosonic degrees of freedom. We therefore 
have to carefully choose observables which are not sen¬ 
sitive to such huctuations. (Their effect on the phase 
diagram of 2cQCD was investigated in Ref. [38].) Exam¬ 
ples of such observables, discussed in this paper, include: 

+ Thermodynamical observables at low temperature 
and high density, where the physics is dominated 
by a Fermi sea of quarks as well as the chiral and 
diquark condensates. Indeed, it is well-known that 
the mean-held BCS pairing theory works quantita¬ 
tively quite well at weak coupling and zero temper¬ 
ature. 

-I- Expectation value of the Polyakov loop from low 
to high temperatures. The Polyakov loop cou¬ 
ples only indirectly to, and therefore is very mildly 
affected by, the colorless diquark degrees of free¬ 
dom [24]. The dominant matter contribution to 
the Polyakov loop thermodynamics comes from the 
colored quarks, which we do take into account. 

Observables for which the diquark huctuations play an 
essential role, and which we therefore deliberately avoid 
discussing in this paper, include: 

— Critical temperature for baryon superhuidity. 
This is the main reason why we prefer not to plot 
the phase diagram of 2cQCD. 

— Baryon number density in the conhned phase with¬ 
out a diquark condensate (that is, at low T as well 
as /tb). 

There are several open questions that we would like 
to understand better and we leave them for future work. 
Firstly, we have not provided any microscopic derivation 
of the chiral twist. Both interaction terms in Eq. (6) are 
expected to arise as we integrate out high-momentum 
modes in 2cQCD to arrive at a low-energy effective de¬ 
scription. In particular, it would be desirable to have an 
explanation for the fact that A < 1. We simply assumed 
this since it leads to the desired phenomenology. 

Secondly, we have not touched upon the fact that in 
the lattice simulations of Ref. [3], a second onset is ob¬ 
served at high ^B, where the thermodynamic quantities 
rise significantly above their SB limits. At present, we 
do not have any explanation for this effect. 


16 


Finally, we wish to make the case that current and 
future lattice data for dense two-color QCD should be 
used to improve our understanding of the real world. To 
this end, we need a dictionary to translate between mod¬ 
els for two-color and three-color QCD. What we have in 
mind is a framework valid for both two and three colors, 
augmented by a mapping of the model parameters. For 
example, one can use the expected scaling based on large- 
Nc arguments. This was done here and in Ref. [24] to fix 
the parameters of the NJL part of the model. The gauge 
sector can be treated in similar manner. Once we have 
a specific analytical model that reproduces or at least 
fits the data for we can extract an improved /re- 

dependent Polyakov-loop potential for QCD simply by 
assuming that its quark-induced part is suppressed by 
a factor of l/N^ relative to the dominant, gluon con¬ 
tribution. The resulting model will allow one to study, 
for example, the interplay of chiral symmetry restoration 
and deconfinement in cold and dense quark matter. 


Appendix A: Chiral model with a Wilson term 


In this appendix, we consider the effects of adding a 
Wilson term to the NJL Lagrangian. 


1. Three-momentum cutoff 

Let us revisit the model of Eq. (6). In the mean- 
field approximation, the Lagrangian describes a pair of 
fermion species with squared masses 

Ml = ^(l + 2 kM ± Vl -b 4 kM) , (AI) 

where M = mo + a. In the limit k —>■ 0, one of the masses 
approaches M, as expected. The other, however, diverges 
as 1/k. In lattice field theory, this mode is interpreted 


as the Wilson doubler which decouples in the continuum 
limit. 

We shall now argue that a naive three-momentum cut¬ 
off A is not a suitable regulator in a theory with a heavy 
doubler. To see this, consider the thermodynamic po¬ 
tential in the vacuum. The thermodynamic potential is 
given by 


f^NJL = ^-4IVcE 

e—± 


d^k 


(A2) 


where eek = \/k‘^ + M|. The gap equation for the chi¬ 
ral condensate is obtained by taking the derivative with 
respect to a. This yields 


a = 


4GA, 

K 



d^k J_ 
(27r)3 tek 



VI + 4kM 


, (A3) 


generalizing Eq. (11). The core of the problem is visi¬ 
ble now. In the limit k —>■ 0, the integrand involving the 
heavy doubler approaches a nonzero constant. Hence the 
heavy fermion gives a nonzero contribution to the gap 
equation even in the limit when its mass goes to infin¬ 
ity. The heavy doubler does not decouple. This can be 
traced to the fact that the Wilson term modifies the spec¬ 
trum of the theory, yet the momentum cutoff is sensitive 
to momenta only, not to masses. We need to use some 
regularization scheme which affects masses directly. 


2. Pauli-Villars regularization 

Implementing the Pauli-Villars (PV) scheme may be 
subtle since covariant regulators in general tend to intro¬ 
duce unphysical artifacts at finite temperature and chem¬ 
ical potential. We therefore first derive the naive expres¬ 
sion for the thermodynamic potential with the Wilson 
term, and subsequently introduce the regulator. This 
will ensure thermodynamic consistency. 

Returning to Eq. (8) and switching to Euclidean space, 
the mean-field thermodynamic potential reads 


f^PNJL = 


gauge 


($) + ^ + ^ _ 2T V 
'' ’ An A\n ^ 


d^k 

(^ 


logdet^ ^(a;„,fe), 


(A4) 


where the determinant of the inverse propagator takes the form 

det^-i = ^[[M - - K'^][{M - KKlf - Kl] + 2\/:^\^[{M - kK'^){M - kKD - ■ K{\ -b |A|4|^ (A5) 


The arrows indicate momenta of red quarks and green 
antiquarks, ^ = (Ar°^,A;), with 

KI = iuJn - 


Furthermore, uJn = (2n -I- l)7rT are the fermionic Mat- 
subara frequencies and a is the constant background 
color gauge field, related to the Polyakov loop variable 
via $ = cos{f3a). It is not easy to factorize the de- 


= vjJn -b /i — ia, 


fi — ia. (A6) 













17 


terminant (A5) into simple factors corresponding to free 
fermionic quasiparticles for general A. In the following, 
we therefore set A = 0 and make sure that /iB is low 
enough so that we do not enter the baryon superfluid 
phase. Note that as a consequence A becomes irrelevant 
for the thermodynamics. 

PV regularization is introduced by making the replace¬ 
ment in Eq. (A4) 

logdet^”^ ^ logdet^”^, (A7) 

3 

where the inverse propagator is dehned analogously 
to Eq. (A5) by 


violet divergences. Two simple choices, referred to as 
the PV2 and PV3 schemes, correspond to j = 0,1,2 
with Cj = {1,—2,1} and to j = 0,1,2,3 with cj = 
{1,-3, 3,-1} [6]. The PV2 scheme removes all the di¬ 
vergences, except a residual logarithmic divergence in the 
thermodynamic potential. This divergence can be elimi¬ 
nated by subtracting the potential at a conveniently cho¬ 
sen reference point. All derivatives of the thermodynamic 
potential such as the gap equation are rendered finite. 
The PV3 scheme, on the other hand, renders the ther¬ 
modynamic potential finite without further subtractions. 
We choose the PV3 scheme for convenience since the PV2 
scheme turns out to be continuous but non-analytic in the 
limit K —>■ 0. 


det 5 -^ = { [{M -KK^f+ - K^] 

2 (A8) 

The coefficients Cj are chosen in order to cancel ultra- 


Since at A = 0 the Dirac determinant (A5) trivially 
factorizes, it is easy to carry out the Matsubara sum, 
leading to a generalization of Eq. (29), valid in the normal 
phase. 


qA=o _ Q 
''‘PNJL “ ^‘gauge V 




e=± ± 


d^k 


E< 

3=0 


Eii +Tlog{l + ) 


(A9) 


where 






M" = — 
2^2 


^1 -b 2kM + e ,\/1 -b 4 kM — 4}k2A2^ 


(AlO) 


The latter generalizes Eq. (Al) for the Wilson masses to 
the PV regulator modes. Before we work out the physical 
consequences, it is useful to pause and explain why we 
introduced the PV regularization via Eq. (A8). Here is 
our consistency check list: 

• The regularization manifestly cancels all diver¬ 
gences. This is easy to see by expanding the loga¬ 
rithm of the regulated determinant (A8) in powers 
of A. Every factor of jA^ is suppressed by 1/A'^ 
for K > 0, or at least l/AT^ for /c = 0. Since the 
thermodynamic potential has a quartic divergence, 
three subtractions are sufficient to make it finite. 

• Eor K = 0, our prescription reduces to the simple 

replacement —>■ M'^+jA^, corresponding to the 

usual implementation of the PV scheme in Lorentz- 
invariant theories. However, this naive replacement 
does not work for nonzero n. 

• Our thermodynamic potential represents a gas of 
quasiparticles with masses Aljf, at finite chemical 
potential. Since the regularization does not affect 
the shift of by /i, it automatically has the “Sil¬ 
ver Blaze” property [36], namely that the physics 
is completely independent of /tb at T = 0 below 


the onset of diquark condensation. Note that this 
property is not necessarily preserved is some PV- 
like regularizations where the whole Gram matrix 
of the Dirac operator is regulated [39]. 

3. Vacuum physics 

As the first step, we have to refit the model parameters 
owing to the different regularization scheme. The same 
values for the chiral condensate, pion decay constant, and 
pion mass given in Eq. (18), yield the parameter values 
in the PV3 scheme 

G = 5.11 GeV■^ 

A = 1129 MeV, (fitted parameters) (All) 
mo = 5.4 MeV. 

The fact that the same physical observables require a 
relatively high cutoff and give a comparably small con¬ 
stituent quark mass in covariant regularization schemes 
is well known [6]. The Wilson coupling n is treated as a 
free parameter. 
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kA 


FIG. 12. Chiral condensate in the vacuum as a function of 
the dimensionless combination kA; all the other parameters 
are fixed by Eq. (All). Note the logarithmic scale of the 
horizontal axis. 

Figure 12 shows the effect of the k coupling on the 
chiral condensate in the vacuum. As expected, the chi¬ 
ral condensate initially grows with k. However, around 
kA « 10“^, the growth stops and further increasing k 
leads to a suppression of the condensate. This is con¬ 
nected to the fact that the PV regulator masses (AlO) 
become complex, and we therefore do not attach physi¬ 


cal significance to this threshold. Large values of n are 
not physical anyway since the masses of the two fermion 
species converge to each other. 

In a similar fashion, one can use Eq. (A9) to study 
the consequences of the Wilson term for thermodynam¬ 
ics quantities. However, we do not perform a detailed 
analysis here, but just remark that it does not affect the 
qualitative conclusions made in this paper. At a quanti¬ 
tative level, we expect k to play a role whenever chiral 
symmetry is important. In particular, the presence of 
an extra heavy fermion species should manifest itself in 
modified thermodynamics at high temperatures. 
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